
import geopandas as gpd
from shapely.geometry import Point, Polygon, MultiPoint
from matplotlib_scalebar.scalebar import ScaleBar
import matplotlib.pyplot as plt

import folium
from folium.plugins import MarkerCluster
import geopandas as gpd

import numpy as np

import gurobipy as gp
from gurobipy import GRB


# In this cell, let's read data

# Data source: granted by SB Fire District GIS Manager
stations = gpd.read_file("Data/Fire_Stations_SBC.shp")

# Data source: https://catalog.data.gov/dataset/tiger-line-shapefile-2019-county-santa-barbara-county-ca-topological-faces-polygons-with-all-ge
county = gpd.read_file("Data/tl_2019_06083_faces.shp") # check out Data folder


county


# Unfortunately, county data contains the ocean part, which we don't want..
    # check it out in QGIS or ArcGIS if you want!

# Let's remove the water part by using the column "LWFLAG". 
    # When LWFLAG column value is L, it means the geometry is part of land,
    # While LWFLAG column value is W, it means the geometry is waterbody
county = county.loc[county['LWFLAG'] == 'L'].reset_index(drop=True)


# Let's unary_union every geometry in county GeoDataFrame
    # Because it contains all the small blocks, but we want the whole county boundary.
county_shape = county.unary_union

    # Note: unary_union returns a GEOMETRY only!


# Plot the two datasets in one plot
# Does it look good..?
A = stations.plot(color="red")

# To plot the GEOMETRY, county_shape with another gpd.GeoDataFrame, stations, 
# let's convert it into gpd.GeoSeries
gpd.GeoSeries(county_shape).plot(ax=A)


# When the plot doesn't look plausible, first thing to do is to look into the CRS.

<Axes: >


# Let's check stations dataset's crs
    # It's in a Projected Coordiante system, where the unit is "foot"
stations.crs

<Projected CRS: EPSG:2229>
Name: NAD83 / California zone 5 (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - California - counties Kern; Los Angeles; San Bernardino; San Luis Obispo; Santa Barbara; Ventura.
- bounds: (-121.42, 32.76, -114.12, 35.81)
Coordinate Operation:
- name: SPCS83 California zone 5 (US survey foot)
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich


# Let's check county dataset's crs
    # It's in a Geographic Coordinate system, where the unit is "degree"
county.crs

<Geographic 2D CRS: EPSG:4269>
Name: NAD83
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. Puerto Rico. United States (USA) - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Hawaii; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands. British Virgin Islands.
- bounds: (167.65, 14.92, -40.73, 86.45)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich


# Let's convert the county dataset to stations crs. 
# Because stations have a planar coordinate system (foot linear unit).

# Actually, what this line does is to create a new GeoDataFrame
    # with one geometry, unary_union of water county blocks, in stations.crs.
county_shape = gpd.GeoDataFrame(geometry=[county.to_crs(stations.crs).unary_union], crs=stations.crs)
county_shape


# Now, let's plot them together.

A = stations.plot(color="red", zorder=100, markersize=10)
county_shape.plot(ax=A, facecolor='None', edgecolor="grey")

# Review for SCALEBAR!
scalebar = ScaleBar(dx=1, scale_formatter=lambda data, unit: f'{data},000 ft', length_fraction=0.25)
A.add_artist(scalebar)

<matplotlib_scalebar.scalebar.ScaleBar at 0x12f821d50>


# How about interactive mapping?

import folium
from folium.plugins import MousePosition
import geopandas as gpd


county_shape_4326 = county_shape.to_crs(4326)

# Calculate the centroid of the polygon
centroid = county_shape_4326.geometry.centroid.values[0].coords[0]

# Create a Folium map centered at the centroid of the polygon
m = folium.Map(location=[centroid[1], centroid[0]], zoom_start=9)

# Add the polygon to the map
folium.GeoJson(county_shape_4326, style_function=lambda x: { 'fillColor': 'transparent'}).add_to(m)

folium.GeoJson(stations.to_crs(4326)).add_to(m)

# Add MousePosition control to display cursor coordinates
#### Use the cursor coordinates to define your Area of Interest! ####
MousePosition().add_to(m)

# Display the map
m

/var/folders/f9/08glrw657c7c1tj4s2fnqw9m0000gn/T/ipykernel_65821/3721739741.py:11: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  centroid = county_shape_4326.geometry.centroid.values[0].coords[0]


# Following 4 float values are given for lab question 1.
# For lab question 2, you need to change it to your own boundary coordinates
max_x, min_x = -119.55734, -119.90547
max_y, min_y = 34.46071, 34.39000

# Create the rectangle geometry, which forms are study region boundary.

boundary = Polygon([(min_x, min_y), (min_x, max_y), (max_x, max_y), (max_x, min_y), (min_x, min_y)])


# Add the boundary to the map
folium.GeoJson(boundary.__geo_interface__, style_function=lambda x: {'color': 'red', 'fillColor': 'transparent'}).add_to(m)

m


# Do we agree with the study region boundary?
# Let's check the area

# let's have it as a GeoSeries, with crs defined (4326)
# And convert it back to feet crs (stations.crs)
boundary = gpd.GeoSeries(boundary, crs=4326).to_crs(stations.crs)[0]

# And retrieve the area value
boundary.area / 2.788e+7
# According to google, let's divide the value with 2.788e+7 to convert square ft to square miles
# The boundary is about 97 square miles

96.89657718898926


# To take a subset of fire stations within our study region, 
    # Let's get the boolean (True/False) array indicating whether each station is within our boundary.

stations.within(boundary).sum()

19


# Remember? We can use .loc[] to take the subset!
    # Then we use reset_index(drop=True) to make the row numbers consecutive.

study_stations = stations.loc[stations.within(boundary)].reset_index(drop=True)

# Check out how many stations we have in our study region.
len(study_stations)

19


# Also, let's take an intersection of county polygon with our study region.


boundary_county = county_shape.intersection(boundary) # the output will be assigned to boundary_county variable
boundary_county

0    POLYGON ((6091967.348 1976096.323, 6091892.338...
dtype: geometry


# Let's define how much fire station can cover!
# I would try 1 mile firstly
COVERAGE_DISTANCE = 1 * 5280 # because 1*5280 feet is 1 mile 

# Draw buffer from each fire station in study_stations GeoDataFrame 
    # Buffer distance is COVERAGE_DISTANCE we've defined
stations_coverage = study_stations.buffer(COVERAGE_DISTANCE)

# Let's plot the county polgyon within our study region, and
    # stations_coverage buffers
A = stations_coverage.plot(facecolor="None")
boundary_county.plot(ax=A, zorder=-1)
plt.axis("off")

# How much area is covered?
    # To get to know, firstly take unary_union of the total station_coverage
total_stations_coverage = stations_coverage.unary_union
    # Then take the intersection with boundary_county, to exclude water area
total_stations_coverage = total_stations_coverage.intersection(boundary_county)
    # Then divide the area of the updated total_station_coverage with the boundary_county.area
coverage = total_stations_coverage.area / boundary_county.area
    # Multiply 100 to get the % value
coverage = coverage * 100

# Check out how much percentage of study region is covered by current fire stations.
    # Do they agree with each other? the coverage % and plot?
print(f'{round(coverage.values[0],2)} % covered')

60.52 % covered


# This is a short example, how to create 2 pair combinations in python.

people = ["Jiwon", "Zhongqi", "Rosemary", "Rita"]

for i in range(len(people)-1):
    for j in range(i+1, len(people)):
        print((i, j), people[i],",", people[j])
        
# Just be familiarized with this syntax :)
    # The combination result makes sense, right?

(0, 1) Jiwon , Zhongqi
(0, 2) Jiwon , Rosemary
(0, 3) Jiwon , Rita
(1, 2) Zhongqi , Rosemary
(1, 3) Zhongqi , Rita
(2, 3) Rosemary , Rita


# How much redundant coverage?

redundant_coverage = 0
    # This two for-loops is the most efficient way to go through all 2 pair combinations of stations_coverage
for i in range(len(stations_coverage)-1):
    # Take one station coverage buffer
    cover_i = stations_coverage[i]
    
    for j in range(i+1, len(stations_coverage)):
        # Take another station coverage buffer which hasn't paired with i
        cover_j = stations_coverage[j]
        print("\t combination: ", i, j)
        
        # With the two stations combination, calculate the intersection area.
            # and sum up!
        redundant_coverage += cover_i.intersection(cover_j).area

print() # To add blank line
print("Redundant Coverage: ", redundant_coverage / 2.788e+7, " square miles")

	 combination:  0 1
	 combination:  0 2
	 combination:  0 3
	 combination:  0 4
	 combination:  0 5
	 combination:  0 6
	 combination:  0 7
	 combination:  0 8
	 combination:  0 9
	 combination:  0 10
	 combination:  0 11
	 combination:  0 12
	 combination:  0 13
	 combination:  0 14
	 combination:  0 15
	 combination:  0 16
	 combination:  0 17
	 combination:  0 18
	 combination:  1 2
	 combination:  1 3
	 combination:  1 4
	 combination:  1 5
	 combination:  1 6
	 combination:  1 7
	 combination:  1 8
	 combination:  1 9
	 combination:  1 10
	 combination:  1 11
	 combination:  1 12
	 combination:  1 13
	 combination:  1 14
	 combination:  1 15
	 combination:  1 16
	 combination:  1 17
	 combination:  1 18
	 combination:  2 3
	 combination:  2 4
	 combination:  2 5
	 combination:  2 6
	 combination:  2 7
	 combination:  2 8
	 combination:  2 9
	 combination:  2 10
	 combination:  2 11
	 combination:  2 12
	 combination:  2 13
	 combination:  2 14
	 combination:  2 15
	 combination:  2 16
	 combination:  2 17
	 combination:  2 18
	 combination:  3 4
	 combination:  3 5
	 combination:  3 6
	 combination:  3 7
	 combination:  3 8
	 combination:  3 9
	 combination:  3 10
	 combination:  3 11
	 combination:  3 12
	 combination:  3 13
	 combination:  3 14
	 combination:  3 15
	 combination:  3 16
	 combination:  3 17
	 combination:  3 18
	 combination:  4 5
	 combination:  4 6
	 combination:  4 7
	 combination:  4 8
	 combination:  4 9
	 combination:  4 10
	 combination:  4 11
	 combination:  4 12
	 combination:  4 13
	 combination:  4 14
	 combination:  4 15
	 combination:  4 16
	 combination:  4 17
	 combination:  4 18
	 combination:  5 6
	 combination:  5 7
	 combination:  5 8
	 combination:  5 9
	 combination:  5 10
	 combination:  5 11
	 combination:  5 12
	 combination:  5 13
	 combination:  5 14
	 combination:  5 15
	 combination:  5 16
	 combination:  5 17
	 combination:  5 18
	 combination:  6 7
	 combination:  6 8
	 combination:  6 9
	 combination:  6 10
	 combination:  6 11
	 combination:  6 12
	 combination:  6 13
	 combination:  6 14
	 combination:  6 15
	 combination:  6 16
	 combination:  6 17
	 combination:  6 18
	 combination:  7 8
	 combination:  7 9
	 combination:  7 10
	 combination:  7 11
	 combination:  7 12
	 combination:  7 13
	 combination:  7 14
	 combination:  7 15
	 combination:  7 16
	 combination:  7 17
	 combination:  7 18
	 combination:  8 9
	 combination:  8 10
	 combination:  8 11
	 combination:  8 12
	 combination:  8 13
	 combination:  8 14
	 combination:  8 15
	 combination:  8 16
	 combination:  8 17
	 combination:  8 18
	 combination:  9 10
	 combination:  9 11
	 combination:  9 12
	 combination:  9 13
	 combination:  9 14
	 combination:  9 15
	 combination:  9 16
	 combination:  9 17
	 combination:  9 18
	 combination:  10 11
	 combination:  10 12
	 combination:  10 13
	 combination:  10 14
	 combination:  10 15
	 combination:  10 16
	 combination:  10 17
	 combination:  10 18
	 combination:  11 12
	 combination:  11 13
	 combination:  11 14
	 combination:  11 15
	 combination:  11 16
	 combination:  11 17
	 combination:  11 18
	 combination:  12 13
	 combination:  12 14
	 combination:  12 15
	 combination:  12 16
	 combination:  12 17
	 combination:  12 18
	 combination:  13 14
	 combination:  13 15
	 combination:  13 16
	 combination:  13 17
	 combination:  13 18
	 combination:  14 15
	 combination:  14 16
	 combination:  14 17
	 combination:  14 18
	 combination:  15 16
	 combination:  15 17
	 combination:  15 18
	 combination:  16 17
	 combination:  16 18
	 combination:  17 18

Redundant Coverage:  19.585007656426257  square miles


# Let's create a function to calculate a redundant coverage!

def calculate_redundant_coverage(coverage_buffers):
    redundant_coverage = 0
        # This two for-loops is the most efficient way to go through all 2 pair combinations of stations_coverage
    for i in range(len(coverage_buffers)-1):
        # Take one station coverage buffer
        cover_i = coverage_buffers[i]

        for j in range(i+1, len(coverage_buffers)):
            # Take another station coverage buffer which hasn't paired with i
            cover_j = coverage_buffers[j]

            # With the two stations combination, calculate the intersection area.
                # and sum up!
            redundant_coverage += cover_i.intersection(cover_j).area

    print() # To add blank line
    print("Redundant Coverage: ", redundant_coverage / 2.788e+7, " square miles")
    return redundant_coverage / 2.788e+7


current_rc = calculate_redundant_coverage(stations_coverage)

Redundant Coverage:  19.585007656426257  square miles


# Check this GeoDataFrame attribute, bounds!
    # It's convenient to find minx, miny, maxx, maxy
boundary.bounds

(5986841.110329453, 1967950.0435882716, 6092251.258445264, 1995487.2553522934)


# Let's save them as variables.

minx, miny, maxx, maxy = boundary.bounds
print(minx, miny, maxx, maxy)

5986841.110329453 1967950.0435882716 6092251.258445264 1995487.2553522934


# Here, we are going to create fishnet points.
    # Fishnet points are points spaced with an equal interval.
    # You'll see when you see the output!

# Define the interval between points 
#### Note: Adjust interval value as needed for lab question 2
interval = 5280/5  # default: 5280/5 feet (1 mile / 5 = 0.2 mile) 

# Create arrays of x and y coordinates using np.arange
x_coords = np.arange(minx, maxx, interval)
y_coords = np.arange(miny, maxy, interval)

# Create a list to store the points
fishnet_points = []

# Generate points for the fishnet
for y in y_coords:
    for x in x_coords:
        fishnet_points.append((x, y))

# Print the number of points generated
print("Number of points in the fishnet:", len(fishnet_points))

fishnet_points = gpd.GeoSeries([Point(pt_cd) for pt_cd in fishnet_points])
fishnet_points.plot(figsize=(15,5))

Number of points in the fishnet: 2700

<Axes: >


### Note: Here, we want the fishnet points only if 
    # they are within our study region, and within the existing stations coverage
    
fishnet_points = fishnet_points.loc[fishnet_points.intersects(boundary_county[0])]
fishnet_points = fishnet_points.loc[fishnet_points.intersects(total_stations_coverage[0])].reset_index(drop=True)

fishnet_points.plot(figsize=(15,5))

<Axes: >


# Overlay the fishnet_points with existing fire stations

A = fishnet_points.plot(figsize=(15,5), color="black", markersize=10)
study_stations.plot(ax=A, color="blue", marker="*", markersize=400)

<Axes: >


# Can we cover the same amount with Fewer stations?
# Let's solve this optimization problem, LSCP !!

num_demands = len(fishnet_points) 
num_facilities = len(study_stations)

D_ij = [[np.sqrt((fishnet_points.x[i] - study_stations.geometry.x[j])**2 + (fishnet_points.y[i] - study_stations.geometry.y[j])**2) 
         for j in range(num_facilities)] for i in range(num_demands)]

N_i = dict(zip(range(num_demands), [[] for _ in range(num_demands)]))

for i in range(num_demands):
    for j in range(num_facilities):
        if D_ij[i][j] <= COVERAGE_DISTANCE:
            N_i[i].append(j)

# Create a new model
model = gp.Model("LSCP")

# Decision variables
# x[i] = 1 if facility i is selected, 0 otherwise
x = model.addVars(num_facilities, vtype=GRB.BINARY, name="x")

# Objective function: minimize the number of selected facilities
model.setObjective(gp.quicksum(x[j] for j in range(num_facilities)), GRB.MINIMIZE)

# Constraints
# Each demand point must be covered by at least one selected facility
for i in range(num_demands):
    model.addConstr(gp.quicksum(x[j] for j in N_i[i]) >= 1, name=f"cover_demand_point_{j}")

# Optimize the model
model.optimize()

Set parameter TokenServer to value "gurobi.geog.ucsb.edu"
Gurobi Optimizer version 11.0.1 build v11.0.1rc0 (win64 - Windows 10.0 (19045.2))

CPU model: 12th Gen Intel(R) Core(TM) i7-12700, instruction set [SSE2|AVX|AVX2]
Thread count: 12 physical cores, 20 logical processors, using up to 20 threads

Optimize a model with 1012 rows, 19 columns and 1391 nonzeros
Model fingerprint: 0x539068d2
Variable types: 0 continuous, 19 integer (19 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective 18.0000000
Presolve removed 1012 rows and 19 columns
Presolve time: 0.01s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.02 seconds (0.00 work units)
Thread count was 1 (of 20 available processors)

Solution count 1: 18 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.800000000000e+01, best bound 1.800000000000e+01, gap 0.0000%


selected_bool = [x[i].x > 0.5 for i in range(num_facilities)]
selected_facilities = study_stations.loc[selected_bool].reset_index(drop=True)
unselected_facilities = study_stations.loc[~np.array(selected_bool)].reset_index(drop=True)

print("Objective value: ", model.objVal, sum(selected_bool))
print("Current running stations: ", len(study_stations))

Objective value:  18.0 18
Current running stations:  19


A = gpd.GeoSeries(MultiPoint(fishnet_points)).plot(figsize=(15,5), color="grey", markersize=5)
unselected_facilities.plot(ax=A, color="blue", legend=True, label="Unselected", marker="*", markersize=300, zorder=10)
selected_facilities.plot(ax=A, color="red", marker="*", markersize=300, legend=True, label="Selected")
selected_facilities.buffer(COVERAGE_DISTANCE).plot(ax=A, edgecolor="red", facecolor="None", linestyle="dashed")
plt.legend()
plt.axis("off")

(5987839.3707838515, 6089675.605640798, 1966194.1857715468, 1997983.11008415)


# Convert GeoPandas GeoSeries to GeoDataFrame
fishnet_gdf = gpd.GeoDataFrame(geometry = fishnet_points, crs = stations.crs).to_crs(4326)

# Create a base map
m = folium.Map(location=[fishnet_gdf.unary_union.centroid.y, fishnet_gdf.unary_union.centroid.x], zoom_start=11)



# Plot fishnet
for idx, row in fishnet_gdf.iterrows():
    folium.CircleMarker(location=[row.geometry.y, row.geometry.x],
                        radius=1,
                        color='grey', fillcolor="grey").add_to(m)

# Plot study stations
for idx, row in unselected_facilities.to_crs(4326).iterrows():
    folium.Marker(location=[row.geometry.y, row.geometry.x],
                   icon=folium.Icon(color='blue'),
                   popup=row['Label']).add_to(m)

# Plot selected facilities
for idx, row in selected_facilities.to_crs(4326).iterrows():
    folium.Marker(location=[row.geometry.y, row.geometry.x],
                   icon=folium.Icon(color='red', icon='star'),
                   popup=row['Label']).add_to(m)

# Plot buffer around selected facilities
for idx, row in selected_facilities.to_crs(4326).iterrows():
    folium.vector_layers.Circle(location=[row.geometry.y, row.geometry.x],
                                 radius=COVERAGE_DISTANCE/3.281, # feet to meter conversion
                                 color='red',
                                 fill=False,
                                 dash_array='10').add_to(m)

# Display the map
m


# Calculate Redundant Coverage of selected facilities
    # Hint: you can use the function we defined!
    
optimized_rc = calculate_redundant_coverage(selected_facilities.buffer(COVERAGE_DISTANCE)) 
    
# Is redundant coverage decreased with optimization?

Redundant Coverage:  14.188941465159218  square miles


len(selected_facilities)

18


# Let's create the spatial efficienty

# Minimum (ideal) number of facilities needed to cover the demands covered by current facilities
# divided by the number of current facilities (fire stations).

spatial_efficiency = len(selected_facilities) / len(study_stations) * 100

print("Spatial Efficienty is ", round(spatial_efficiency, 2), " %")

# Higher efficienty, closer to the ideal number
# Lower efficienty, more redundancy, further from the ideal number

Spatial Efficienty is  94.74  %


print("Coverage Distance (mi): " , COVERAGE_DISTANCE) # Double Check this value before putting into table ** 
print("Current Stations #: " ,len(study_stations))
print("Optimized Stations #: ", len(selected_facilities))
print("Current Redundant Coverage (sqmi): ", current_rc)
print("Optimized Redundant Coverage (sqmi): ", optimized_rc)
print("Spatial Efficiency (%): ", spatial_efficiency)

Coverage Distance (mi):  5280
Current Stations #:  19
Optimized Stations #:  18
Current Redundant Coverage (sqmi):  19.585007656426257
Optimized Redundant Coverage (sqmi):  14.188941465159218
Spatial Efficiency (%):  94.73684210526315

Lab 2: Fire Station Service Coverage Analysis¶

Overview¶

Objectives¶

Tools and Technologies¶

Dataset¶

Lab Structure¶

Step 1. Data Exploration and Preprocessing¶

Step 2. Analysis of Existing Facility Placement¶

Step 3. Optimization of Facility Placement¶

Location Set Covering Problem¶

Step 4. Evaluation of Spatial Efficiency Metrics¶

Final Report¶

	TFID	STATEFP10	COUNTYFP10	TRACTCE10	BLKGRPCE10	BLOCKCE10	SUFFIX1CE	ZCTA5CE10	UACE10	PUMACE10	...	METDIVFP	CNECTAFP	NECTAFP	NCTADVFP	LWFLAG	OFFSET	ATOTAL	INTPTLAT	INTPTLON	geometry
0	257492722	06	083	002211	2	2004	None	93454	None	08301	...	None	None	None	None	L	N	220098	+34.9497690	-120.3695124	POLYGON ((-120.37608 34.95466, -120.36430 34.9...
1	257492275	06	083	002211	2	2003	None	93454	None	08301	...	None	None	None	None	P	N	5308	+34.9546966	-120.3766078	POLYGON ((-120.37705 34.95521, -120.37608 34.9...
2	265426787	06	083	002211	2	2005	None	93454	None	08301	...	None	None	None	None	P	N	5945	+34.9550536	-120.3774668	POLYGON ((-120.37794 34.95483, -120.37769 34.9...
3	257492266	06	083	002211	2	2005	None	93454	None	08301	...	None	None	None	None	P	N	149887	+34.9585978	-120.3839902	POLYGON ((-120.39081 34.96311, -120.37769 34.9...
4	219296658	06	083	003102	1	1258	None	93436	None	08302	...	None	None	None	None	L	N	255	+34.4938573	-120.4754525	POLYGON ((-120.47580 34.49391, -120.47569 34.4...
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
19184	219313538	06	083	002906	2	2004	None	93117	79282	08303	...	None	None	None	None	L	N	104502	+34.4500980	-119.8355449	POLYGON ((-119.83790 34.44967, -119.83785 34.4...
19185	219312610	06	083	002005	2	2001	None	93455	79417	08301	...	None	None	None	None	L	N	115	+34.8613367	-120.4128687	POLYGON ((-120.41296 34.86140, -120.41282 34.8...
19186	219314638	06	083	002809	1	1050	None	93436	None	08302	...	None	None	None	None	P	N	2023	+34.6601685	-120.2962177	POLYGON ((-120.29654 34.66031, -120.29643 34.6...
19187	219312971	06	083	002502	3	3045	None	93434	None	08302	...	None	None	None	None	P	N	16783	+34.9595483	-120.6316773	POLYGON ((-120.63511 34.95954, -120.63416 34.9...
19188	219313149	06	083	001800	1	1416	None	93252	None	08302	...	None	None	None	None	P	N	410	+34.8472318	-119.4805966	POLYGON ((-119.48072 34.84728, -119.48065 34.8...